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Commensurate-incommensurate transitions are ubiquitous in physics and are often accompanied 
by intriguing phenomena. In few-layer graphene (FLG) systems, commensurability between honey- 
comb lattices on adjacent layers is regulated by their relative orientation angle 6, which is in turn 
dependent on sample preparation procedures. Because incommensurability suppresses inter-layer 
hybridization, it is often claimed that graphene layers can be electrically isolated by a relative twist, 
even though they are vertically separated by a fraction of a nanometer. We present a theory of 
interlayer transport in FLG systems which reveals a richer picture in which the specific conduc- 
tance depends sensitively on 0, single-layer Bloch state lifetime, in-plane magnetic field, and bias 
voltage. We find that linear and differential conductances are generally large and negative near 
commensurate values of 9, and small and positive otherwise. 



Experimental advances in the fabrication of graphene- 
based structures [l|,[l] have now provided researchers with 
a multitude of systems that have strikingly distinct elec- 
tronic properties. By engineering the substrate un- 
derlying exfoliated samples identifying exfoliated 
fragments with folds @, or controlling epitaxial growth 
conditions 0, ttie size and shape of the honeycomb 
lattice arrays^^ and the number of graphene lay- 
ers and their orientations can all be varied. This struc- 
tural diver sity nourishes hopes for a future carbon-based 
electronics I llj with band-structure and transport char- 
acteristics that can be tailored for different types of ap- 
plications. 

FLG has advantages over single-layer-graphene be- 
cause it has a larger current-carrying capacity and be- 
cause its electronic properties are sensitive to more en- 
gineer able system parameters In nature it appears 
in a variety of stacking arrangements, the most com- 
mon being Bernal and rhombohedral sequences which can 
form three dimensional lattices. It has been understood 
for some timejTsI that in graphite can depart from 
Bernal values. With some interesting exceptions [sl. 
most recent studies of inter-layer twists in FLG have fo- 
cused on samples grown on SiC[l^. In particular Hass 
et. al. have demonstrated that orientational disorder 
is normally present in carbon- face SiC epitaxial FLG 
samples [ig]. The present work is motivated primarily 
by the need to achieve a more complete understanding of 
transport in these graphitic nanostructures, which cur- 
rently appear to provide the most promising platform 
for applications. 

In a bilayer system, the relative rotation angle 9 can be 
classified as either commensurate or incommensurate [T7| . 
In the former case the misaligned bilayer system still 
forms a crystal, albeit one with larger lattice vectors and 
more than four atoms per unit cell. Commensurability 
occurs at a countably infinite set of orientations; but the 
probability that a randomly selected orientation angle is 
commensurate vanishes. The energy bands of commensu- 
rate twisted multilayers disperse approximately linearly 
with momentum |18l420l | , except at energies very close to 
the Dirac point. However, the Dirac velocity is reduced 
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FIG. 1: Interlayer (RC) equilibration rate as a function of 
twist angle 6. These results were calculated for two layers 
with equal carrier densities (n — 5 x lO^^cm"^) and EfT = 
3, where ef is the Fermi energy and r is the isolated-layer 
Bloch state lifetime. The relaxation rate is dominated by 
separate features that appear near every commensurate angle, 
but differ in strength by many orders of magnitude. The tails 
of individual features have been cut-off in this plot in order to 
reveal weaker features that will emerge in more ideal bilayers. 
Except near 6 = 0, the equilibration rate is surprisingly slow 
for two layers separated by an atomic length scale. 



compared to that of a single layer swtem especially for 
rotation angles close to 0° or 60° [13 ■ The linear 
Dirac-like dispersion contrasts with the approximately 
quadratic dispersion found in a Bernal stacked bilayer 
systemim. Incommensurate bilayers are not crystalline 
and therefore their electronic properties cannot be ana- 
lyzed using Bloch's theorem. 

Here we develop a theory of the vertical transport 
properties of twisted FLG samples which is valid in the 
incoherent transport limit [1^. We show that the spe- 
cific linear conductance between misaligned layers is en- 
hanced over a small but finite range of twist angles near 
those that produce relatively short period commensurate 
structures, that the conductance peak angles shift with 
in-plane magnetic field -B|| , and that the peaks become 
narrower and stronger when the isolated layer Bloch state 
lifetime r increases. The differential conductivity tends 
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to be negative near commensurate conductance peaks 
and positive otherwise. Typical theoretical results for 
the dependence of the interlayer equilibration rate on 6 
are presented in Fig. [T] In the following we first explain 
the analysis which supports these statements and then 
discuss some implications for FLG electronics. 

Studies of transport between weakly coupled two- 
dimensional (2D) electron systems have a long historyfU, 
[25I in semiconductor hetero junctions systems. In that 
case epitaxial tunnel barriers are responsible for nearly 
perfect 2D momentum conservation, which then helps to 
make vertical transport a powerful probe of electronic 
properties. Our theory of vertical transport in FLG is 
similar to the successful semiconductor heterojunction 
theory [23|. We derive an expression for tunneling current 
/ vs. bias voltage V by using a 7r-orbital tight-binding 
model, approximating inter-layer hopping processes at 
leading order in perturbation theory, and accounting for 
the inevitable presence of a finite disorder potential which 
limits the life-times of Bloch states in each layer. These 
steps lead to 
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where gs = 2 accounts for spin degeneracy, AiQ(k, w) is 
the spectral function for band a and layer i, Upi is the 
Fermi distribution function for layer i, and T^pi is the 
tunneling matrix element between isolated layer Bloch 
states with band and crystal momentum labels, |ka) and 
|p'/3). The sums over k and p' may be taken over the 
unrotated and rotated Brillouin zones respectively. We 
derive Eq. ([T]) in section 2 of the Supplementary Infor- 
mation , where we justify its neglect of disorder vertex- 
corrections. In our calculations, A is approximated by 
a Lorentzian function with full-width-half-maximum h/r 
centered on the band energy cta (k) . (Hereafter h = 1 and 
length is measured in units of Uc = 1.42A, the carbon- 
carbon distance in graphene.) Eq. ([T]) is valid in the 
weak tunneling regime in which T is smaller than life- 
time broadening 1/r, allowing coherent tunneling pro- 
cesses to be neglected. This condition is satisfied in typ- 
ical samples except at rotation angles very close to 0° or 
60°. 

In a twisted bilayer system the tunneling matrix ele- 
ment depends strongly on the relative orientation of the 
two graphene sheets. The honeycomb lattice vectors of 
the rotated layer R' are related to those of the unrotated 
layer R by R' = M{9)Ti + d. Here M is the trans- 
formation matrix for rotations in the lattice plane and 
d is a translation vector. Corresponding rotations oc- 
cur in reciprocal space so that eia(jp) — (p') when 
p' = M{9)p. Commensurability is determined only by 
M, but linear translations of one layer relative to the 
other do modify T, and hence the tunneling current. 

The magnitude of T depends on the 7r-orbital interlayer 
hopping amplitudes of our tight-binding model which we 



estimate using a simple two center approximation scheme 
explained in section 1 of the Supplementary Information. 
We find that 
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where fio is the area of a unit cell. Here Gi and G2 are 
summed over reciprocal lattice vectors, primed wavevec- 
tors are rotated, s and s label the two triangular hon- 
eycomb sublattices centered at positions Tg, and a^J is 
the sublattice projection of the jka) Bloch state in the 
unrotated layer. In Eq.([2|), which is derived in section 1 
of the Supplementary Information, is the 2D Fourier 
transform of the finite-range inter-layer hopping ampli- 
tude. As we will explain, the interlayer conductance and 
the layer equilibration rate are proportional to |tkp val- 
ues for |k|'s that are larger than the Brillouin-zone scale 
(except for « 0°, 60°). Because the inter-layer distance 
is already larger than the carbon-carbon distance within 
a layer, these |ikP values tend to be both extremely small 
and extraordinarily sensitive to details of the inter-layer 
tunneling model that are otherwise inconsequential. 
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FIG. 2: Fermi circles in an extended zone scheme. The 
blue (large) and red (small) circles correspond to the Fermi cir- 
cles in the unrotated and rotated layers respectively. The area 
enclosed by the circles is proportional to the carrier density. 
Conductance contributions occur when the Fermi circles inter- 
sect and are much larger when the intersection occurs closer 
to the origin of momentum space. The Brillouin-zone bound- 
ary connects the centers of the inner shell of blue circles as 
indicated by the dashed lines in the 6 = 17° panel. 

We have used Eqs.dl]) and © to evaluate interlayer 
currents as a function of rotation angle 9, carrier density, 
bias voltage, and disorder strength. Since for typical elec- 
tronic densities the temperature T is much less than the 
Fermi temperature we focus on T = hereafter. It is 
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helpful to first focus on the linear conductance 

= ^El^kp'l'^i"(k>^FM2/5(p',e.). (3) 
^ kp' 

The equilibration rate plotted in Fig[T] was obtained by 
viewing the bilayer as a leaky capacitor and ignoring any 
screening by graphene a orbitals. This model yields an 
RC circuit with time constant Trc related to the conduc- 
tance by Gj A — 0.027 t^^ where A is the layer area in 
rr? , G is measured in Siemens, and r^c in seconds. Apart 
from a change in scale, Fig [T] can then be viewed as a plot 
of the interlayer conductance. We find that the tunneling 
conductance increases abruptly near commensurate an- 
gles, that the height of the peaks scales linearly with epT 
(for epT > 1), and that the peaks narrow as r increases. 
The discontinuous jumps of log(G') in Fig (T] are artificial 
and result from a numerical procedure in which momenta 
k and p' in Eq.® are restricted to the vicinity of the 
Fermi energy. This procedure suppresses the tails of all 
commensurate features, allowing more minor features to 
be revealed. In practice the conduction tails correspond- 
ing to highly commensurate structures will dominate G 
over a range of angles that depends on t. Limited by 
computational power we considered « 75 meV in 
FigH] however in epitaxial graphene the lifetime can be 
more than an order of magnitude longer [2^. An accu- 
rate theory of the conduction-peak tails would require 
a reliable theory of the isolated-layer spectral function 
tails. 

Why is the tunneling conductance enhanced at com- 
mensurate rotation angles? To understand the relation 
between interlayer current and commensurability it is il- 
luminating to plot the Fermi surfaces of both layers, pe- 
riodically extended in momentum space by adding re- 
ciprocal lattice vectors to the crystal momenta of the 
electrons. As we see in Eq.([2]), allowed interlayer tun- 
neling processes are diagonal in this generalized momen- 
tum. The left panels in Figl5] corresponds to the in- 
commensurate rotation angles B — 17°, 26° whereas the 
right panels correspond to the commensurate angles near 
Q = 21.8°, 27.8°. We use different Fermi surfaces sizes for 
clarity; similar considerations apply independent of the 
sign or magnitude of the carrier density ratio. The key 
feature to notice in these plots is that at commensurate 
rotation angles some Fermi spheres overlap. Overlaps 
of circles centered on the extended Dirac points, always 
accompany commensurate real-space structures because 
the set of extended Dirac points forms a momentum space 
honeycomb lattice that differs from the real space honey- 
comb lattice only by a scale factor and by a rotation. If 
overlaps occur in real-space, they also occur in momen- 
tum space. Notice that this property holds only when 
the Brillouin-zone corners are extended to fill momentum 
space; if the Dirac point occurred elsewhere in the iso- 
lated layer Brillouin-zone, the dependence of inter-layer 
conductance on 9 would be quite different. 

The overlap of extended Dirac points does not fully ex- 
plain the conductance peaks at finite density, since Fermi 



energy states at finite carrier density are displaced from 
the Dirac point. The nesting between Fermi surfaces al- 
luded to in Fig 12] actually depends not only on commen- 
surability, but also on the fact that for typical carrier 
densities the Fermi surface is well approximated by a 
circle centered on the Brillouin-zone corners. For equal 
densities then, matching Dirac points implies complete 
Fermi surface nesting (see Fig 13]). When the two-layers 
have different densities, the peak conductance will not 
occur at the nesting angle; instead the conductance will 
have a double-peak structure with features offset to both 
sides of the commensurate angle. 

Commensurate rotation angles can be classified as ei- 
ther inter-valley or intra-valley. In the former the two 
Dirac points fco and k'^ that coincide in the extended 
momentum picture are associated with different valleys 
(in the aligned bilayer) whereas for intra-valley rotation 
angles they belong to the same valley. An inter-valley 
commensurate rotation is illustrated in Fig 13] 

Away from commensurate angles the energy differ- 
ence between states which have the same extended mo- 
mentum is typically much larger than the Fermi energy, 
and the spectral function width 1/r (see left panels of 
Figl2]). The conductance is therefore very small away 
from the commensurate-angle peaks. The Dirac-like lin- 
ear spectrum of an ideal commensurately twisted bilayer 
does not, as is commonly stated, indicate that the ideal 
twisted layers are decoupled. At commensurate angles 
the perfect crystal wavefunctions near the Dirac point are 
in fact coherent equal weight contributions from the two 
layers. In the limit of large in-plane Bloch state lifetimes, 
the conductance becomes very large and eventually the 
incoherent transport picture will fail. 




FIG. 3: Nesting of Dirac cones at commensurability. For 
commensurate rotation angles every momenta state on the 
rotated Fermi circle is mapped onto a momenta state of an 
unrotated Fermi circle. 

As wc have explained, vertical transport at commen- 
surability is dominated by processes in which an electron 
tunnels from a momentum near a Dirac point of one layer, 
to a momentum that is the same distance from a Dirac 
point of the other layer. Since carrier densities per unit 
cell are always small, we can replace ik+G in Eq.(15]) by 
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^ko+G where is the Dirac point momentum. We then 
find that the conductance peak can be expressed as the 
product of geometry-related and phase space factors: 

G«i?(0„d)^Ai(k,ep)Ai(k,ep), (4) 

k 

where 9c is the commensurate orientation. R{9ctA) de- 
pends mainly on the value of i|kD+G| at which the ex- 
tended Dirac points overlap (see FigH]), while the re- 
maining phase space factor is identical to the one that 
appears in the theory of coupled quantum wells [23|. For 
equal densities in the two layers, the Fermi surfaces nest 
precisely. For pure rotations R can be calculated analyt- 
ically. For inter-valley commensurate rotation angles we 
find that 

G(^c)=-4,.,.^^^fi^. (5) 

Here Eg is the energy difference between the top con- 
duction band and bottom valence band of the twisted 
bilayer at the Dirac point, and we assumed that e-pT > 1. 
In section 4 of the Supplementary Information we derive 
Eg. (1301) and obtain a similar formula for intra- valley ro- 
tation angles. In addition we numerically verify that the 
conductance changes only by a factor of order unity as d 
is varied across the unit cell. Eq.(|30|) therefore provides a 
good estimate for G regardless of the relative translation 
between the two layers (see section 3 of the Supplemen- 
tary Information). 

When the densities differ, Fermi circles in different lay- 
ers begin to overlap near 9 = 6c only after a momentum- 
space relative shift Q equal in magnitude to the differ- 
ence of the two Fermi wavevectors. As in semiconductor 
double-wellsfH, m 113, a shift Q = z x ed^/l^, where 
III is the magnetic length, can be accomplished in a bi- 
layer with layer separation d± by applying an in-plane 
magnetic field B^^e. For graphene bilayers, however, a 
relative momentum space shift can also be achieved by 
rotation, as is clear from Figl^l For small departures 
from commensurability Q Ri {6 — 9c) z x (kn + G). For 
equal densities, both rotations and in-plane fields dra- 
matically suppress the conductance peak when vQ > l/r 
where v is the band velocity of graphene. For example, 
for n = 4 • lO^^cm"^ and r = 50 fsec [3l|, the conduc- 
tance peak should nearly completely disappear at 0.15 
Tesla. FLG should therefore provide a palette on which 
gate voltages, in-plane magnetic fields, and rotations can 
be mixed to produce a rainbow of interrelated and ex- 
traordinarily strong magnetic-field and strain sensitive 
resistance effects. 

In general the interlayer conductance G{9) is peaked 
whenever any extended Fermi surface overlap occurs at 
reasonably small reciprocal lattice vectors. The degree of 
overlap can be parameterized by q* , the minimum sepa- 
ration between extended Dirac points of the rotated and 
unrotated layers for a reciprocal lattice vector truncation 
chosen to refiect the scale on which t^, falls off. For a clean 
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FIG. 4: The minimum separation between extended Dirac 
points q* as a function of rotation angle 6. 



system, tunneling conductance at equal density is appre- 
ciable as long as q* k, \0 — ^dlk-f- G| < 2ki,. Because 
2kp in FLG electronic systems is always small compared 
to reciprocal lattice vector scales, the conductance peaks 
are invariably sharp when plotted as a function oi 9. As 
an example q* « 6.39|6' — 9c\ in the vicinity of 9c = 27.8° 
for the reciprocal lattice vector illustrated in Figl2j In 
FigHl g*, minimized over the first two G-shells, is plot- 
ted as a function of angle. Overlap between the Fermi 
spheres of the two layers will therefore persist over the 
angle range for which q* is smaller than ^kp- 

Electronic structure calculations for ideal commensu- 
rate bilayers demonstrate that Eg decreases very rapidly 
as the number of atoms per unit cell increases [17|. Eg = 
780 meV for a unit cell of 4 atoms, and already less 
that 1 meV for a unit cell of 100 atoms. It is there- 
fore plausible that conductance tails that correspond to 
the few most lowest order commensurate angles (e.g. 
9c = 0°, 21.8°, 27.8°, 32.2°, 38.2°, 60°) will dominate G at 
every rotation angle. Eqs.Q and ([50]) should therefore 
be interpreted as a lower bound for the conductance at 
higher order commensurate 0c's. 

We now turn to the non- linear I ~V oi twisted bilayer 
graphene. At zero temperature 

I{S,V) = ^ VlTfepf r du:A^{k,u)A2(p',u:+eV). 

(6) 

We numerically find that the I-V curves at commensurate 
and incommensurate angles are drastically different. At 
relatively small bias voltage the currents corresponding 
to commensurate angles are several orders of magnitude 
larger than their incommensurate counterparts. On the 
other hand negative differential conductances invariably 
appear at commensurate angles, whereas dl/dV tends to 
be small and positive at incommensurate angles. 

In classic tunneling experiments, a bias voltage induces 
an equal electric potential difference between the layers. 
Total energy conservation then implies kinetic energy 
changes equal to eV upon tunneling. Since, as we have 
explained, the allowed tunneling processes at commen- 
surate angles are between states with the same kinetic 
energy bias voltages tend to decrease tunneling currents. 
Following the same approximations that led to Eq.(|31) we 
can capture this effect mathematically by expressing the 



5 



interlayer current in product form: 

/oo 
duj[n^{w)^n^{u: + V)] 

CP 

X Y.A^{k,u:)Ap{k,uj + eV) (7) 

K 

In Eq.( [T]) we have allowed for both intraband and in- 
terband tunneling at large biases. As long as eV < 
tunneling between conduction bands dominate / when 
both layers are n-type. In this intermediate non-linear 
regime the two Lorentzian shaped spectral functions in 
Eq.© overlap only weakly and 

/(0.,d)«G(e„d)^^^^. (8) 

for CpT > 1. Negative differential conductance occurs 
when eVr > 1. For incommensurate twist angles, crys- 
tal momenta conservation can not be sustained at the 
Fermi surface. Increasing V unblocks processes in which 
tunneling occurs between states with different kinetic en- 
ergies and leads to a slow increase of the tunneling cur- 
rent with a complex dependence on tq and H/t. For 
eV > ep, the current increases monotonically with V for 
both commensurate and incommensurate twist angles. 
The commensurate tunneling current has a sharp rise at 
eV = 2tp due to momentum conserving processes allowed 
at high bias voltage in which a valence band electron in 
one layer tunnels to the conduction band of the oppo- 
site layer. For commensurate angles it follows from ([7]) 
that these inter-band processes eventually dominate the 
tunneling current and that 

^R-^Q{V -2e^)V (9) 

to leading order in 1/Vt. Here 9 is the Heaviside step 
function. The finite temperature corrections to Eas. (|8l9p 
are exponentially small in T/cp. 

The extension of our theory to FLG is straightforward 
in the linear regime. In the simplest case each layer is ro- 
tated with respect to its neighbors sufficiently to drive 
the system into an incoherent transport regime. The 
weak links between layers then act like classical resistors 
which appear in series in vertical transport. The resis- 
tance of each link depends on the rotation angle between 



layers and on the densities in both layers. We anticipate a 
very rich and complex behavior of FLG in the non-linear 
regime. The negative differential conductivities are likely 
to give rise to steady state multistability and to chaotic 
temporal response, as occurs in semiconductor multiple- 
quantum- well systems [s^l- A more complicated scenario 
could arise in turbostratic graphene. There the entire 
layered structure is composed of a set of coherent multi- 
layer substructures, characterized by either a Bernal or 
an AA stacking sequence. Weak links which play a dom- 
inant role in limiting vertical conductance appear due to 
occasional twists. The calculations for the resistance of 
each twisted interface closely follow those outlined above 
for the two-layer case when supplemented by a band in- 
dex for the various 2D energy bands of a coherent sub- 
structure. 

One application of our theory is to assess whether or 
not twisted graphene layers are effectively isolated from 
an electrical point of view. The equilibration time be- 
tween layers that are spatially uniform but out of equi- 
librium is plotted in Fig. [1] and is very long compared 
to characteristic electronic time scales for rotation an- 
gles far from important commensurabilities, near 10° for 
example. The steady-state equilibration length between 
separately contacted layers can be estimated by equating 
inter-layer conductances, which are proportional to sam- 
ple area, with the intra-layer conductance per square. 
For the commensurate angle 9c = 21.8°, for example, the 
sample area at which they are identical is approximately 
0.04/im^. As evident from FiglU the corresponding areas 
for small rotation angles near the AA and AB stacking 
sequences are even smaller. For small rotation angles, 
the two layers are therefore strongly coupled. 

Finally we remark that the extraordinary sensitivity of 
the tunneling conductance to the twist angle found here 
suggests that misaligned graphene bilayers might be use- 
ful as ultra-sensitive strain gauges or pressure sensors [ssj 
which are widely used in biological, mechanical and op- 
tical systems. 

The authors acknowledge support from CERA, SWAN 
and the Welch Foundation and helpful conversations with 
W. de Heer, R. Duine, P. First, D. Goldhaber-Gordon, 
R. Lifshitz, and E. Tutuc. 



Supplementary Information 

I. THE TUNNELING MATRIX ELEMENTS 

The interlayer hopping terms in a 7r-band tight-binding Hamiltonian for twisted graphene bilayers depend in general 
on the positions of all carbon atoms. Our analysis of inter-layer conductance and equilibration is based on a simple 
two-center model in which the interlayer hopping parameter between two sites, t[r), depends only on the planar 
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projection of their separation r. In the main text we used an equation, derived below, which relates the inter-layer 
hopping amplitudes of twisted bilayers to tq, the two-dimensional Fourier transform of t(r). 

One strategy which can be used to estimate t(r) is to assume functional forms for the distance dependence of the 
Slater-Koster tpp^r and tpp^^ hopping functions [341] . and then fit them to accurately known parameters of untwisted 
bilayers. We have explored this approach, following the procedures adopted in Refs. [13, [111, but have concluded that 
it tends to underestimate hopping amplitudes near the Dirac points of twisted bilayers. We have therefore decided to 
obtain numerical estimates by directly fitting an ansatz for tq to obtain 

iq = e-"(«''-)\ (10) 

where io = 2 eVA^, a — 0.13, 7 = 1.25, and d± — 3.34A is the distance between the layers. The value used for to is 
the average of values implied by the models in Refs. [l3, [11]. Since to is the sum of all inter-layer hopping parameters, 
it should be estimated reliably by any parameterization that uses accurate values for the largest hopping parameters. 
We fix a and 7 so that the values of the ideal bilayer gaps are accurate for the lowest order commensurate structures. 
These are proportional to tk^ {& = 0° and 9 = 60°) and ^6.4/0, {0 = 21.8° and 9 = 38.2°) where = 1.42A is 
the carbon-carbon distance in single layer graphene. See details in Sec. IIVI below. We fit the energy gaps to values 
extracted from the ab initio calculations by Shallcross et a/.[l3]- Note that these values of tq characterize short- 
distance roughness in the inter-layer hopping landscape which survives Fourier transformation at large wavevectors, 
which is not simply related to typical inter-layer hopping strengths. The energy gaps that we obtain at 9 — 21.8° 
using the real space parameterizations of tppa{r) and tppT^{r) in Refs. [13, [111 are both substantially smaller than the 
ah initio gaps of Shallcross et al. [Tt} . 

We now derive the expression for the hopping amplitude between Bloch states in twisted bilayers that is used in 
the main text. The Bloch state in layer j with crystal momentum k and band index a can be written as 



where A and B label the two triangular honeycomb sublattices. 



^ la 
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71 



(12) 



and 0k is the phase of the inter-sublattice hopping term in the single-layer tight-binding model. For nearest neighbor 
hopping within the planes 8k = arg e^^'^^^ where the 5j are the three vectors connecting an atom with its nearest 
neighbors. The Bloch state projection on sublattice s is 

l^W) = 4=^eW-)lr, + R), (13) 



^ R 



where jr, -I- R) is a site-representation basis function of the tight-binding model. In Eq.([T3|) R is a triangular lattice 
vector, N is the number of unit cells in the system, and we choose Ta — Q and Tb equal to the vector connecting the 
two atoms within a unit cell. 

The relative orientation of the two layers can be described by a rotation matrix M{9) and a translation vector d. 
Therefore every Bloch wave function in the second layer is related to a Bloch wave function in the first layer by 

= (14) 

with ]R -I- Ts) in layer 1 replaced by ]R' + r^) in layer 2, r' = Mr + d for all positions and k' — A/k. Using primes to 
indicate layer 2 variables and invoking the two-center approximation for the inter-layer tunneling amplitude, 

{Ts + -R\H,„ter\T' s + R') = t(Ts + R - t', - R'), (15) 



we find that 

ss' R,iR.2 



Expression (2) in the main text is obtained by Fourier expanding ^(r) and summing over the lattice vectors. 
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II. VERTEX CORRECTIONS 



The general expression for the tunneling current 

X /mG?^„ (kN , ko , u:)ImGl^, {p'a,p'^,u: + eV) (16) 

is obtained using second order perturbation theory [sB]. In Eq. (fT6| rij is the Fermi distribution in layer j, G^^^ is the 
retarded Green function in layer j that correspond to the propagation of a charge carrier from band a to band 7, the 
rotation angle is 9, and V is the bias voltage. The over-line denotes disorder averaging. As in the main text, primed 
variables are associated with the rotated layer. Since disorder breaks translation invariance, the Green functions are 
not diagonal is the momentum representation. When the disorder averages can be performed independently for the 
two-layers, translational invariance is recovered and Eq. (ITB)) reduces to Eq.(l) of the main text. 

We average over disorder using the self-consistent Born approximation in which correlations between the layers 
appear as a vertex-correction ladder diagram sum (see Fig|5]). For simplicity we assume white noise disorder and 
characterize the correlation between the disorder potentials in the two layers by 7 = ni{UiU2) where rii is the 
concentration of impurities and Uj is the disorder potential in layer j. For aligned bilayers with short range tunneling 
we find that 

where /3 = ni{U'j). As evident from Eq.(fT7|) the tunneling conductance diverges if the disorder potentials of the 
two layers are perfectly correlated. These strong correlations are likely in a graphene bilayer because of the small 
distance between the layers. The divergence of G indicates the breakdown of perturbation theory, i.e. it invalidates 
the incoherent theory we use in this work. A similar scenario arises for tunneling between coupled semiconductor 
quantum wells[23] when their disorder potentials are strongly correlated. 

We now show that vertex corrections are important only at very small values of the rotation angle 9. The physical 
origin of this behavior is twofold. First, the relevant correlation in the twisted case is between the disorder potential 
in one layer and a spatially rotated counterpart in the other layer. For any finite range disorder correlation length, 
these two disorder potentials are independent making 7 in Ea. (|17p considerably smaller. Second, the divergence in 
the conductance appears due to tunneling between identical states. However, for incommensurate angles the wave 
vectors of the initial and final states in a tunneling process substantially differ making (3 in Ea. ([T7|) considerably 
larger. In the following paragraphs we explain how this latter behavior is captured in a diagrammatic perturbation 
theory description of a disordered system. 

We first focus on the tunneling conductance for aligned layers [9 — Q). At zero temperature 

G(0) = ^ E T^kop^^kt;;,^"^^?, JkN, ko, 6,)/mG5,,(p;,, pi^, £,). (18) 
The conservation of crystal momentum in expression (2) for T^^, implies that po — kg and that pn = kN. For 




FIG. 5: Self consistent Born approximation. A bubble diagram with ladders. 

CpT > 1 interband transitions are inhibited so that a = 7 and [3 = 5. Due to the spinor form of the wave functions 
each disorder line contributes [1 -I- cos{9k — 6'fc )]7/2 to the ladder diagram. To evaluate n("), the ladder diagram 
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with n disorder lines, we first integrate over the angular variables using 

/ ^ cos(0fc, - Oq) cos{dk, ~Og)^l cos{0k, ~0k,)- (19) 
Then using J^(0) — 2'Kv-pT where 

HQ)^Y.^Ul.^)GtM + Q,Lo) (20) 

we integrate over the radial direction. In obtaining J^(0) we have replaced the energy dependent density of states by 
VY^ its value at the Fermi energy. We find that for n>\ 

n(") = G5'(fco)G^(fco) 

where /? = l/nv-pT and /z, = i?, A. For n>2 the Green functions in one layer are retarded and those of the other 
layer are advanced. We now sum XI*-") to infinite order in n. While the sum can clearly be carried for a general 
tunneling matrix element T^^ the basic physical idea is more transparent for short range tunneling. Therefore in the 
calculations below we assume T^k = t \s momentum independent in which case we recover Eq. ()17p . 

We now address the role played by vertex corrections is twisted bilayers. As in the main text our discussion excludes 
the vicinity of 6* = 0°, 60° for which i > l/r. The procedure outlined above for calculating G can be repeated for any 
rotation angle 0. For a rotated bilayer it follows from Eq.(2) in the main text that kg — pj, = G'^ — Gi = Q where 
Q ~ Ikn + Gi||6' — 6c\- Expression (flT)) can then be used for a rotated bilayer as well if fi is replaced by 

p^^pF{Q)/F{Q). (22) 

Because is a monotonically decreasing function of its argument and because Q is comparable in size to the Dirac 
momentum /3q ^ j5. 




FIG. 6: Dependence of conductance per unit cell G on translation d for 9 = 27.8°, e-pT = 3, and n = 5 ■ lO^'^ cm 



III. DEPENDENCE OF TUNNELING CURRENT ON TRANSLATION 

Commensur ability depends only on the relative rotation of the two graphene layers. Nevertheless linear translation 
of one layer with respect to the other will change the tunneling current. In Figl6]the conductance at = 27.8° is 
plotted as a function of d for a bilayer with n = 5 • 10"'^^ cm~^ in each layer and CpT = 3. The dependence of the 
tunneling current on d is captured by the phase factor exp[— i(k + Gi) • d] ~ exp[— i(kD + Gi) • d] in the tunneling 
matrix element Tkp'. When summed over k the result is a rapid spatial variation on the lattice constant scale. 



1 + cos(6lfeo 



^G5'(fcjv)G^(A:w) 



(21) 
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illustrated in Fig|6]due to the exp[— i(kD + Gi) • d] factor, modulated by a slower variation on the Fermi wavelength 
scale. 

As illustrated in Fig.l the conductance peaks appear symmetrically around 6 ~ 30°, but the height of a peak with 
< 30° does not necessarily equal the height of the corresponding peak at 0' — 60° ~ 9. In fact, the relative height of 
the two peaks depends on d. An AA stacking sequence can be transformed to Bernal stacking either by a pure rotation 
with 9 = 60° or by a translation with d = (1, 0). Since the latter transformation docs not influence commensurability 
any commensurate angle 9c of the AA stacked bilayer is a commensurate angle of the Bernal stacked bilayer. The 
conductance peaks then lie symmetrically with respect to 9 = 30° since if 9 is commensurate so is its inverse. 



IV. CONDUCTANCE FOR COMMENSURATE ANGLES FOR d = 0. 

For commensurate angles the conductance can be approximated by Eq.(4) in the main text. The integration over 
the overlap of the two spectral functions 

Ai(k,eF)Ai(k,eF) = [2 + 27rei^r + 4eFT arctan (2eFT)] (23) 

^ — ' znv' 

k 

is independent of the rotation angle and the entire dependence of G on the relative alignment of the two layers is in 
Rfj,v{0, d). We now evaluate i? for d = 0. 

At the Dirac point the intra-layer Hamiltonian vanishes and we have contributions only from interlayer tunneling: 

Ho=(^t[) (24) 

where each element is a 2 x 2 block for the two 7r-bands in each layer. Using a representation of sublattice sites in 
each layer we find that 

(2\a\ \ /2|a| 0\ . . 

Here and T° correspond respectively to intra-valley (S=same) and inter-valley (D=different) rotation angles as 
explained in the main text, and (j) = 0, ±60 depends on 9c: e.g. ^(0°) = 0, (/)(27.8°) = 60° and ^(38.2°) = -60°. If 
the hopping amplitude tq decreases fast enough with momentum so that only the first G-shell significantly contribute 
to the tunneling matrix 

|a| ^ i.5%t£i (26) 

"0 

where fig is the area of a unit cell and Gi is the wavevector which produces the smallest q extended-zone Dirac-cone 
overlap as explained in the text. In our model Eq. ((26|) is satisfied for all commensurate angles except for = 0°, 60° 
for which \a\ = 1.67ikD+Gi/^^o- DiagonaHzing Hq yields — ±2\a\ (both doubly degenerate) and = 0, 0,±2|a|. 
In both cases the energy gap between the top conduction band and bottom valence band is therefore Eg = A\a\. 

To find R we assume that T is well approximated by Eq. (pS]) for finite momentum states in the vicinity of the Dirac 
points. We verified this assumption numerically for low densities. We first focus on inter-valley rotation angles. In 
the eigenstate representation 



It follows from Eqs. (l25l27|) that 



Consequently, 



and 



r,;': = 4-4er,;fv. (27) 

^ g^(2e,+e,)|^|j^ (28) 
/ |^|^^'^(^^.)P = |a|/ (29) 



G°(6le, (i = 0) = Agygs^ J \^i [2 + 27reFT + 4eFr arctan (2eFr)] . (30) 



10 



Expression (5) is obtained in the e-pT > 1 limit. 
Similarly for intra-valley rotation angles 



»:<tr4> 1 (31) 



2 



It then follows that 



f.s cos2 (0 - f ) sin^ v</^ - \ (oo) 



2 y 2 , 

and that the conductance is 

G^(0c, d^0)= Ag.gs^^i^ cos" f - ^) [2 + 2neFr + 46^^ arctan (2ej.T)] . (33) 

Note that inter-band resonant conduction (which occurs when the carrier densities in the two layers are opposite) 
has the same form as its intra-band counterpart for inter-valley rotation angles. For the inter-band conduction at 
inter- valley rotation angles the cos function in Eq. p3p should be replaced by a sin. Interestingly, the ratio 

depends only on the twist angle. For example, AG(27.8°) = 1.94 in accord with the numerical results depicted in 

Figm 

Using the momentum dependent T matrices we can find the bands in the vicinity of the Dirac points. For inter-valley 
rotations we find four non-degenerate bands 



El=±\lel + 2\a\^±2\lel\a\' + \a\\ (35) 
At low energies ^ a 

E^i=±l^ , E^2^±2\a\±^ (36) 
2m* 2m* 

where m* — \a\/v'^. For the intra-layer rotations 



El = ±^el + 4|a|2 ± 4|a|efe cos (0 - ^c/2) (37) 

At low energies ^ a 

El = ±2|a| ± v*k (38) 

where v* = ucos((/) — 9c/2). Deviations from expressions (j25p for T result in trigonal warping in a bilayer system. 
More elaborate studies of the spectrum are needed to determine whether such effects are important in a rotated 
bilayer system as well. 
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